Francesco Pirotti
2025-09-16
CRS - coordinate reference system
https://r.geocompx.org/spatial-class#crs-intro
https://r-spatial.org/book/02-Spaces.html#sec-crs
meuse dataset
https://cloud.r-project.org/web/packages/gstat/vignettes/gstat.pdf
metodi di interpolazione geografica
https://r-spatial.org/book/12-Interpolation.html
Stimare i valori di una variabile nel dominio dello spazio (eventualmente anche del tempo)
Cosa si intende per variabilità spaziale?
Useremo il dataset “meuse” e le librerie che vedete sotto
library(tidyverse)
library(stats)
library(ggplot2)
library(terra)
library(gstat)
library(kableExtra)
library(sp)
data(meuse)
class(meuse)## [1] "data.frame"
kable(head(meuse, 5), caption = "Dataset *meuse*") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))| x | y | cadmium | copper | lead | zinc | elev | dist | om | ffreq | soil | lime | landuse | dist.m |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 181072 | 333611 | 11.7 | 85 | 299 | 1022 | 7.909 | 0.0013580 | 13.6 | 1 | 1 | 1 | Ah | 50 |
| 181025 | 333558 | 8.6 | 81 | 277 | 1141 | 6.983 | 0.0122243 | 14.0 | 1 | 1 | 1 | Ah | 30 |
| 181165 | 333537 | 6.5 | 68 | 199 | 640 | 7.800 | 0.1030290 | 13.0 | 1 | 1 | 1 | Ah | 150 |
| 181298 | 333484 | 2.6 | 81 | 116 | 257 | 7.655 | 0.1900940 | 8.0 | 1 | 2 | 0 | Ga | 270 |
| 181307 | 333330 | 2.8 | 48 | 117 | 269 | 7.480 | 0.2770900 | 8.7 | 1 | 2 | 0 | Ah | 380 |
Vogliamo un modello per “spiegare” la variabilità e dunque poter stimare (interpolare) al meglio dove non abbiamo informazioni.
Stimare valori di una variabile applicando un modello creato conoscendo valori misurati (e.g. campionamenti). E’ ragionevole adottare un criterio spaziale, ovvero usare un criterio di spazio (il valore del punto sarà probabilmente più simile ad una misura vicina che una misura più lontana, a parità di altre condizioni). Da qui la geostatistica.
Richiamo concetti statistici fondamentali: varianza, covarianza, correlazione, autocorrelazione
Richiamo concetti base su sistemi di riferimento coordinate (CRS), geografiche e proiettate
Fondamenti statistici dell’analisi spaziale dati e variabili regionalizzate
Metodi deterministici di stima per interpolazione spaziale (IWD)
Vediamo come valutare alcuni momenti di una distribuzione - media e varianza ci danno il primo e secondo momento. Z è la variabile aleatoria, nel caso del dataset MEUSE possiamo scegliere la concentrazione di Zinco.
\[M(Z)=\mu_Z=\mathbb E[Z]\]
\[M(Z)=\mu_Z=\frac{1}{n}\sum_{j=1}^n{Z_j}\]
\[Var(Z)=\sigma^2_Z=\mathbb{E}\Big[\big(Z-\mathbb{E}[Z]\big)^2\Big].\]
\[Var(Z)=\sigma^2_Z=\frac{\sum_{j=1}^n(Z_j-M_Z)^2}{n}\]
Formula più pratica
\[Var(Z)=\sigma^2_Z=\mathbb{E}[Z^2]-\mathbb{E}[Z]^2\]
\[Var(Z)=\sigma^2_Z=M(Z^2)-M(Z)^2 \]
Calcoliamo la varianza con tre formule diverse:
MZ<-mean(meuse$zinc)
VarZ <- var(meuse$zinc)
## media quadratica dei quadrati
MZ2<-mean(meuse$zinc^2)
VarZ <- MZ2 - MZ^2
VarZ <- mean((meuse$zinc-MZ)^2)
plot( meuse$y, meuse$zinc, xlab="Coord Y", ylab="Zinco", pch="+")
abline(a = MZ,0 )
rect(min(meuse$y), MZ-VarZ^0.5, max(meuse$y),
MZ+VarZ^0.5, col = "#ff000045", border = "#ffffff00")hist(meuse$zinc)
rug(meuse$zinc)
rect(MZ-VarZ^0.5, 10000,
MZ+VarZ^0.5, 0, col = "#ff000045", border = "#ffffff00")
abline(v=MZ, col="red", lwd=2)La covarianza tra due variabili casuali \(Z_1\) e \(Z_2\) misura quanto esse tendono a variare insieme. Se positiva variano entrambe nella stessa direzione, se negativa in direzioni opposte.
\[\mathrm{Cov}(Z1,Z2)=\mathbb{E}\Big[\big(Z1-\mathbb{E}[Z1]\big)(Z2-\mathbb{E}[Z2]\big)\Big]\]
La covarianza dipende dalle unità di misura delle variabili → non è standardizzata. Per questo si usa spesso la correlazione, che normalizza la covarianza tra -1 e 1.
La correlazione = 1 equivale ad una covarianza che ha valore uguale alla varianza.
## [1] TRUE
Analizziamo un gradiente di temperatura simulato
load("oggetti.rda")
plot(lat_seq, temps, type = "l", col = "blue", lwd = 2,
xlab = "Latitude (°N)", ylab = "Temperature (°C)",
main = "Gradiente Temperatura: Italia → Danimarca")
points(lat_seq, temps, pch = 16, col = rgb(0,0,1,0.3))La covarianza equivale alla varianza se la calcoliamo confrontandola a se stessa, ma se trasliamo di un dato intervallo (“lag”) e la ricalcoliamo, la covarianza non diventa zero, diminuisce gradualmente….
## [1] 6.150969
## [1] 5.122209
## [1] 4.972588
## Provate con una variabile senza correlazione
# rand <- runif(200)
# cov(rand,rand)
# cov(rand[2:200],rand[1:199])
# cov(rand[3:200],rand[1:198])Questo è un concetto molto importante da comprendere.
Più aumentiamo la distanza del lag, come è ragionevole pensare per dati spazialmente correlati, più diminuisce il valore di covarianza (si allontana dal valore della varianza == covarianza ).
Se si abbassa per poi rialzarsi, può indicare un pattern ricorrente
(provate sostituendo i valori di temps con un valore ciclico come
temps <- sin((1:200)/10) .
nm <- length(temps)
plot(temps[1:200], temps[1:200], pch="x", cex=0.5, xlab="Temperatura", ylab="Temperatura" )
points(temps[1:199], temps[2:200], pch="x", cex=0.5, col="red" )
points(temps[1:198], temps[3:200], pch="x", cex=0.5, col="green" )Non ci crediamo? Rifacciamolo per un dato casuale con media e varianza uguale a quella del nostro dato:
t.rand <- rnorm(nm, mean = MZ, sd = VarZ^0.5)
plot(t.rand[1:(nm-1)], t.rand[2:nm], pch="x", cex=0.5, xlab="km", ylab="Temperatura" )
cc=1
covs.rand<-list("0"=cov(t.rand, t.rand), "1"=cov(t.rand[1:(nm-1)], t.rand[2:nm] ) )
for(lag in seq(6, 36, 10)){
covs.rand[[as.character(lag)]] <- cov(t.rand[1:(nm-lag)], t.rand[(lag+1):nm])
points(t.rand[1:(nm-lag)], t.rand[(lag+1):nm], cex=0.5, pch="x", col=cc )
cc=cc+1
}Relazione lineare tra due variabili, cosiddetta di Pearson. Semplicemente normalizza la covarianza per le deviazioni standard delle due variabili moltiplicate tra loro:
\[-1\le\rho_{Z1Z2}=\frac{Cov(Z1Z2)}{\sqrt{Var(Z1)}\cdot\sqrt{Var(Z2)}}=\frac{\sum_{i=1}^n(Z1-{\mu_{Z1}})(Z2-{\mu_{Z2}})}{\sqrt{\sum_{i=1}^n(y_i-{\mu_{Z1}})^2}\sqrt{\sum_{i=1}^n(y_i-{\mu_{Z2}})^2}}\le+1\]
L’autocorrelazione misura quanto una variabile è correlata con sé stessa quando viene considerato uno spostamento nel dominio dei valori.
In una sola dimensione — sia essa tempo o spazio (ad esempio nel nostro transetto di temperatura SUD → NORD) — questo spostamento viene detto lag, traducibile come sfasamento.
L’analisi dell’autocorrelazione si effettua tramite la funzione CCF (Cross-Correlation Function), che nel caso specifico valuta la correlazione della serie con sé stessa a diversi lag.
acf() in R ti mostra l’autocorrelazione ai vari lag (sfasamenti).
Lag = 1 → confronta ogni punto con il successivo.
Lag più grandi → confrontano punti sempre più distanti lungo il transetto.
#ccf(temps, temps, plot = F)
#ccf(temps, temps)
# Autocorrelazione
acf(temps, main = "Autocorrelazione della temperatura")Rappresentare informazioni associate alla loro posizione nello spazio
Due modelli: il modello vettoriale ed il modello raster.
Ogni valore è georiferito, ovvero associato ad un punto in un sistema di riferimento, quindi delle coordinate.
L’accuratezza può variare a seconda dello strumento usato per la misura.
L’accuratezza dipende dalla scala, che determina la tolleranza agli errori.
I dati vettoriali hanno tre primitive geometriche di base: punti linee e poligoni. La base è sempre il punto con una coordinate, in 2d o 3d.
Ci sono decine di formati digitali per i dati vettoriali.
Per un riferimento vedete la pagina QUI.
Ad oggi gli strumenti in ambiente R di riferimento per leggere dati vettoriali è la libreria sf.
Il formato raster è una semplificazione/discretizzazione dello spazio in celle/pixel/ nodi a passo costante, quindi una “griglia” georiferita mediante uno o più punti. Il passo della griglia è spesso identico nelle due direzioni (ma non necessariamente), e ne determina la risoluzione spaziale.
Può essere 3D (voxel)
Ci sono decine di formati di file che rappresentano dati raster. Per un riferimento vedete la pagina QUI.
Ad oggi la libreria di riferimento per leggere questi dati è la libreria terra, e stars
Attenzione alla gestione della memoria in quanto i dati raster possono diventare molto pesanti dal punto di vista di spazio di memoria (1000 pixel x 1000 pixel è un raster “piccolo” ma già fornisce 1e6 elementi.
Sono delle convenzioni - si stabilisce un sistema con una origine ed un orientamento degli assi tra loro perpendicolari - Coordinate planimetriche o polari - Sistemi locali o globali
Prima di tutto bisogna definire un riferimento fisico - la forma della superficie definita come “interfaccia solida” è “complicata” (è elastica e subisce l’effetto di svariate forze che la deformano)… e neanche stabile nel tempo!
Si prende come convenzione la forma che avrebbe la superficie equipotenziale (livello del mare) se non ci fossero terre emerse - lo si chiama GEOIDE. E’ la superficie di riferimento per comprendere dove la gravità spinge l’acqua. NB è in corso con il recente piano di telerilevamento nazionale il calcolo di un nuovo nuovo geoide ufficiale nazionale (vedi qui)
Il geoide è fondamentale per comprendere fenomeni legati alle quote, ma molto complesso da trattare. Si semplifica con l’ellissoide, una rappresentazione trattabile matematicamente.
Ci sono tanti ellissoidi che sono stati stimati per avvicinarsi il più possibile alla realtà del geoide
Latitudine e longitudine, coordinate polari.
Coordinate planimetriche, su un piano. Ma la terra non è piatta!
Soluzione: accettiamo deformazioni ma le teniamo sotto controllo (scala, accuratezza, tolleranza)
Si proiettano i punti nella “sfera” su un piano.
Vedi documentazione di QGIS QUI
Le coordinate di Prato della Valle: Lat: 45.398557, Long. 11.876460 trasformate nei vari sistemi. TEST - prova a consultare il documento IGM e trasformare la coordinata nel sistema “RDN2008 / Zone 12 (N-E)” e vedere come differisce dal sistema “RDN2008 / Italy zone (N-E)” … sono uguali?
library(sf)
lat<-45.398557
long<-11.876460
# punti in classe SF
pt <- st_sfc(st_point(c(long, lat)), crs = 4326)
# UTM fuso 32
pt_utm <- st_transform(pt, 32632)
#
pt_ecef <- st_transform(pt, 4978)
# le coordinate nei vari sistemi
geographic = sf::st_coordinates(pt)
planimetric = st_coordinates(pt_utm)
geocentric = st_coordinates(pt_ecef)variabili misurate in punti e correlate spazialmente
Possono essere trattate come campioni di un fenomeno continuo (non discreto)
Presentano autocorrelazione spaziale
Vedi slides iniziali
Media, mediana, varianza, deviazione standard
Distribuzione dei dati (istogrammi, boxplot)
Identificazione di outlier - per il trattamento rigoroso delle osservazioni (non lo trattiamo)
Differenze tra dati indipendenti e dati spazialmente correlati
I metodi deterministici non considerano la componente aleatoria (incerta, casuale) dei dati, ma si basano su assunzioni deterministiche per stimare i valori in punti non campionati - ovvero non considerano la struttura spaziale, nessuna autocorrelazione e non forniscono incertezza.
Stimano il valore Z* in un punto dove non abbiamo una misura usando le misure note vicine, ma come?
\[Z* = F(x,y) = \sum_{i=2}^n w_i Z_i\] dove il peso del punto \(1\) dipende dalla distanza \(h\) di un \(i-esimo\) numero di punti \(n\)
Semplicemente quello più vicino
\[ Z*= Z_i \] Dove \(Z_i\) è il vicino più prossimo e \(Z_*\) il punto da interpolare
poligoni di Thiessen danno l’area di influenza.
\[ Z*= Z_i \] dove il peso è indicato da Sibson
\[ Z_i(\mathbf{x})=\frac{A(\mathbf{x}_i)}{A(\mathbf{x})} \]
dove \(A(x)\) è il volume della nuova cella centrata in x e A(xi) è il volume dell’intersezione tra la nuova cella centrata in x e la vecchia cella centrata in xi.
interpolatore esatto: i valori dei dati originali vengono mantenuti nei punti di riferimento.
crea una superficie liscia priva di discontinuità.
è interamente locale, in quanto si basa su un sottoinsieme minimo di posizioni dei dati che esclude le posizioni che, pur essendo vicine, sono più distanti di un’altra posizione in una direzione simile.
è adattivo dal punto di vista spaziale e si adatta automaticamente alle variazioni locali della densità dei dati o della disposizione spaziale.
privo di ipotesi statistiche (puramente deterministico).
può essere applicato a insiemi di dati molto piccoli, poiché non è basato su dati statistici.
privo di parametri, quindi semplice
\[Z* = F(x,y) = \sum_{i=2}^n w_i Z_i\] dove il peso dipende dalla distanza \(h\)
\[w_i = \frac {h_i^{-p}}{\displaystyle \sum_{j=i}^n h_j^{-p}}\]
e da una potenza \(p\) che va inserita dall’utente e può essere oggetto di ottimizzazione.
library(gstat)
data(meuse.grid)
data(meuse)
## molto importanti i passaggi sotto
coordinates(meuse) = ~x+y
coordinates(meuse.grid) = ~x+y
gridded(meuse.grid) = TRUE
###
zinc.idw = idw(zinc~1, meuse, meuse.grid)## [inverse distance weighted interpolation]
## [1] "SpatialPixelsDataFrame"
## attr(,"package")
## [1] "sp"
La libreria terra riconosce l’oggetto “sp” - posso aprirlo con qualsiasi sofware GIS! Nota: non vengono specificati sistemi di riferimento.
Conoscendo l’area, vicino a Maastricht - Paesi Bassi - applichiamo un SR probabile e verifichiamo aprendo in GIS
terra::crs(zinc.idw.terra) <- "EPSG:28992"
# terra::writeRaster(zinc.idw.terra, "zincIDW.tif", overwrite=T)Eseguite la procedura sopra per interpolare la variabile “lead” (piombo). Salvate il GeoTIFF con il sistema di riferimento corretto nella cartella condivisa QUI dentro la cartella “giorno1” nominando il file cognome_nome.tif